In [1]:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline
在线性回归的例子中，我们假设 $y|x;\theta \sim N(\mu, \sigma^2)$；而对于逻辑回归，我们假设 $y|x;\theta \sim b(1, \phi)$，也即伯努利分布。本节我们将看到，这两个模型都是广义线性模型的特殊形式。
本节包括以下内容：
指数分布族定义如下： $$ P(y;\eta) = b(y)exp(\eta^TT(y)-a(\eta)) $$
其中，$\eta$ 称为自然参数 natural parameter，$T(y)$ 称为充分统计量 sufficient statistic（在我们接下来将考虑的分布中，$T(y)=y$），$a(\eta)$ 称为对数配分函数 log partition function。$e^{-a(\eta)}$ 通常作为使分布 $p(y;\eta)$ 之和(对于离散型)或积分（对于连续型）为 $1$ 的正则化常数。
固定函数 $T, a, b$，可以得到以 $\eta$ 为参数的一族分布。随着 $\eta$ 的改变，我们可以得到同族的不同分布。
伯努利分布和正态分布都是指数分布族中的一族分布。
以伯努利分布为例： $$ \begin{split} p(y; \phi) &= \phi^y(1-\phi)^{1-y} \\ &= exp(ylog\phi + (1-y)log(1-\phi)) \\ &= exp((log(\frac{\phi}{1-\phi}))y + log(1-\phi))) \end{split} $$
这里，我们令 $\eta=log(\frac{\phi}{1-\phi})$，可以解出得到$\phi=\frac{1}{1+e^{-\eta}}$，正是sigmoid函数。
此外，令 $T(y) = y, a(\eta) = -log(1-\phi) = log(1+e^\eta), b(y) = 1$，即可得到上式。注意这里 $\eta$ 是常数，因而 $\eta^T = \eta$。
考虑正态分布，在我们推到线性回归的过程中，方差 $\sigma^2$ 的值对最终的 $\theta$ 和 $h_\theta(x)$并没有影响。为简化计算，我们设 $\sigma^2 = 1$： $$ \begin{split} p(y; \mu) &= \frac{1}{\sqrt{2\pi}}exp(-\frac{1}{2}(y-\mu)^2) \\ &= \frac{1}{\sqrt{2\pi}}exp(-\frac{1}{2}y^2) \cdot exp(\mu y-\frac{1}{2}\mu^2) \end{split} $$
令 $\eta=\mu, T(y)=y, a(\eta)=\frac{\mu^2}{2}=\frac{\eta^2}{2}, b(y)=\frac{1}{\sqrt{2\pi}}exp(-\frac{y^2}{2})$，可得上式。
此外，诸如多项式分布、泊松分布、伽马分布、指数分布、Beta分布、狄利克雷分布等等，都属于指数分布族。
给定一组特征 $x$，我们希望预测随机变量 $y$。在这种情况下，构造一个广义线性模型，大致需要以下几个假设：
$y|x;\theta \sim ExponentialFamily(\eta)$。也即是说，给定 $x$ 和 $\theta$，假设$y$ 服从以 $\eta$ 为参数的某一种指数分布族。
给定 $x$，我们的目标是预测 $T(y)$ 的期望值。在大多数例子中，$T(y)=y$。所以这意味着我们希望假设函数 $h$ 的输出符合条件 $h(x) = E[y|x]$。注意到，线性回归和逻辑回归都满足这个条件。
自然参数 $\eta$ 和输入 $x$ 是线性相关的，即 $\eta = \theta^Tx$。（或者，如果 $\eta$ 是向量的话，$\eta_i=\theta_i^Tx$）
其中，第三点某种程度上可以视作我们的一个算法设计选择，而并非本质意义上的假设。通过这三个假设，构造出的广义线性模型，通常易于学习，对于预测不同分布的 $y$ 也常常十分有效。
在二元分类问题中，随机变量 $y \in \{0, 1\}$ 是离散型随机变量，假设其服从伯努利分布 $b(1, \phi)$，则其期望为 $\phi$。于是假设函数为： $$ \begin{split} h_\theta(x) &= E[y|x; \theta] \\ &= \phi \\ &= \frac{1}{1+e^{-\eta}} \\ &= \frac{1}{1+e^{-\theta^Tx}} \end{split} $$
每行的推导过程同最小二乘法。因此，使用sigmoid函数，实际是在伯努利分布假设之上定义广义线性模型的结果。
在多元分类问题中，随机变量 $y \in \{1, 2, ..., k\}$ 依然是离散型随机变量，但现在目标值的范围超过了两个。我们假设其服从多项分布，并根据广义线性模型来建模。首先需要证明，多想分布也属于指数分布族。
多项分布可以用 $k$ 个参数 $\phi_1, ..., \phi_k$ 来描述，每个参数表示取对于结果的概率。但注意到，离散型随机变量取各个可能值的概率和需要等于 $1$，使用 $k-1$ 个参数 $\phi_1, ..., \phi_{k-1}$，同时令 $\phi_k=1-\sum_{i=1}^{k-1}\phi_i$ 就足够了。为了描述方便，我们依然沿用 $\phi_k$ 的写法，但记住这并不是多项分布的一个参数。
我们定义 $T(y) \in \mathbb{R}^{k-1}$ 为：
$$ T(1) = \begin{bmatrix} 1\\ 0\\ 0\\ \vdots\\ 0 \end{bmatrix}, T(2) = \begin{bmatrix} 0\\ 1\\ 0\\ \vdots\\ 0 \end{bmatrix}, T(3) = \begin{bmatrix} 0\\ 0\\ 1\\ \vdots\\ 0 \end{bmatrix}, ..., T(k-1) = \begin{bmatrix} 0\\ 0\\ 0\\ \vdots\\ 1 \end{bmatrix}, T(k) = \begin{bmatrix} 0\\ 0\\ 0\\ \vdots\\ 0 \end{bmatrix} $$和线性回归、逻辑回归都不同，在这里 $T(y) \neq y$。并且 $T(y)$ 在这里是一个 $k-1$ 维的向量，而不再是实数。我们用 $(T(y))_i$ 来表示向量 $T(y)$ 的第 $i$ 个元素。
我们再引入另一个方便的标记符号，函数 $1\{\cdot\}$ 在参数为真时取 $1$，否则取 $0$。这样上面的向量可以简化为 $(T(y))_i = 1\{y=i\}$；并且 $E[(T(y))_i] = P(y=i) = \phi_i$。
按照上面引入的这些符号，我们有： $$ \begin{split} p(y;\phi) &= \phi_1^{1\{y=1\}}\phi_2^{1\{y=2\}}\cdot\cdot\cdot\phi_k^{1\{y=k\}} \\ &= \phi_1^{1\{y=1\}}\phi_2^{1\{y=2\}}\cdot\cdot\cdot\phi_k^{1-\sum_{i=1}^{k-1}1\{y=i\}} \\ &= \phi_1^{(T(y))_1}\phi_2^{(T(y))_2}\cdot\cdot\cdot\phi_k^{1-\sum_{i=1}^{k-1}(T(y))_i} \\ &= exp((T(y))_1log(\phi_1) + (T(y))_2log(\phi_2) + \cdot\cdot\cdot\ + (1-\sum_{i=1}^{k-1}(T(y))_i)log(\phi_k)) \\ &= exp((T(y))_1log(\frac{\phi_1}{\phi_k}) + (T(y))_2log(\frac{\phi_2}{\phi_k}) + \cdot\cdot\cdot\ + (T(y))_{k-1}log(\frac{\phi_{k-1}}{\phi_k})) + log(\phi_k)) \\ &= b(y)exp(\eta^TT(y)-a(\eta)) \end{split} $$
其中： $$ \begin{split} \eta &= \begin{bmatrix} log(\frac{\phi_1}{\phi_k})\\ log(\frac{\phi_2}{\phi_k})\\ log(\frac{\phi_{k-1}}{\phi_k}) \end{bmatrix} \\ a(\eta) &= -log(\phi_k) \\ b(\eta) &= 1 \end{split} $$
这样我们可以确定多项回归属于指数分布族。接下来我们需要求解 $\phi_i$，为方便起见，我们定义 $\eta_k = log(\frac{\phi_k}{\phi_k}) = 0$，这样： $$ \begin{split} e^{\eta_i} &= \frac{\phi_i}{\phi_k} \\ \phi_ke^{\eta_i} &= \phi_i \\ \phi_k\sum_{i=1}^ke^{\eta_i} &= \sum_{i=1}^k\phi_i=1 \\ \phi_k &= \frac{1}{\sum_{i=1}^ke^{\eta_i}} \\ \phi_i &= \frac{e^{\eta_i}}{\sum_{i=1}^ke^{\eta_i}} \end{split} $$
上式的最终结果，就是softmax函数。有了以上的结果，我们可以开始构造广义线性模型： $$ \begin{split} h_\theta(x) &= E(T(y)|x; \theta] \\ &= E \begin{bmatrix} 1\{y=1\}\\ 2\{y=2\}\\ \vdots\\ 1\{y=k-1\}\\ \end{bmatrix} \\ &= \begin{bmatrix} \phi_1\\ \phi_1\\ \vdots\\ \phi_{k-1}\\ \end{bmatrix} \\ &= \begin{bmatrix} \frac{exp(\theta_1^Tx)}{\sum_{j=1}^kexp(\theta_j^Tx)}\\ \frac{exp(\theta_2^Tx)}{\sum_{j=1}^kexp(\theta_j^Tx)}\\ \vdots\\ \frac{exp(\theta_{k-1}^Tx)}{\sum_{j=1}^kexp(\theta_j^Tx)}\\ \end{bmatrix} \end{split} $$
换言之，我们的假设函数将会对所有 $k$ 计算 $p(y=i|x;\theta)$（注意，实际上模型只会计算 $k-1$ 种输出的概率）。可以看到，Softmax回归实际上是逻辑回归的一种泛化形式。按照类似的方法去构造最大似然函数，之后用梯度下降或牛顿法求解 $\theta$ 即可。
In [ ]: